
% This function returns our proposed short and simple confidence intervals
% (SSCIs)

% - Y is the vector of coefficient estimates
% - Sigma is the estimate of the asymptotic variance of Y (Sigma/n)
% - index is the index of the variable of interest
% - type is the "type" of CI, i.e., upper one-sided (1), two-sided (0), or
% lower one-sided (-1)
% - alpha is the nominal level; default is set to 0.05
% - bounded is a vector indicating whether the corresponding entry is
% bounded below by the corresponding entry in bounds (1), unbounded (0), or
% bounded above by the corresponding entry in bounds (-1); default is a
% vector of 1s (i.e., all nuisance parameters are bounded from below)
% - bounds is the vector of bounds; default is a vector of 0s (i.e., all
% nuisance parameters are bounded (from below or above) by zero)

function CI = SSCI(Y,Sigma,index,type,varargin)

%% Assign default values for alpha, bounded, and bounds if they are not provided %%
optargs = {0.05 ones(length(Y),1) zeros(length(Y),1)};

optargs(1:length(varargin)) = varargin;

[alpha, bounded, bounds] = optargs{:};

%% Error handling %%

% Check if Y is a vector
if ~isvector(Y)
    error('Y needs to be a vector!');
    return;
else
    l = length(Y);
end

% Check if Sigma is a matrix
if ~ismatrix(Sigma)
    error('Sigma needs to be a matrix!');
    return;
else
    [nrows,ncols] = size(Sigma);
end

% Check the dimension of Sigma
if nrows ~= l || ncols ~= l
    error('Sigma does not have the right dimension!');
    return;
end

% Check whether Sigma is a variance-covariance matrix
is_var_cov = 1;

if issymmetric(Sigma)
    try
        chol(Sigma);
    catch
        is_var_cov = 0;
    end 
else
    is_var_cov = 0;
end

if is_var_cov == 0
    error('Sigma is not a variance-covariance matrix!');
    return;
end

% Check if index is an integer in the range between 1 and the length of Y
if prod(size(index)) ~= 1
    error('index needs to be scalar!');
    return;
elseif rem(index,1) ~= 0
    error('index needs to be an integer!');
    return;
elseif index < 1 || index > l
    error(['index needs to be between 1 and ' num2str(l) '!']);
    return;
end

% Check if type is in {-1,0,1}
if prod(size(type)) ~= 1
    error('type needs to be scalar!');
    return;
elseif ~ismember(type,[-1 0 1])
    error('type needs to be in the set {-1,0,1}!');
    return;
end

% Check if alpha is in {0.01,0.05,0.1}
if prod(size(alpha)) ~= 1 
    error('alpha needs to be scalar!');
    return;
elseif ~ismember(alpha,[0.01 0.05 0.1])
    error('alpha needs to be in the set {0.01,0.05,0.1}!');
    return;
end

% Check if bounded is a vector of the right length
if ~isvector(bounded)
    error('bounded needs to be vector!');
    return;
elseif length(bounded) ~= l
    error('bounded does not have the right length!');
    return;
end

% Check if bounds is a vector of the right length
if ~isvector(bounds)
    error('bounds needs to be vector!');
    return;
elseif length(bounds) ~= l
    error('bounds does not have the right length!');
    return;
end

%% Read in polynomial coefficients (copy-pasted from paper)
if abs(type) == 1
    if alpha == 0.01
        poly_coef = [2.3476 2.5073 -19.6229 65.0489 -122.0242 112.9814 -40.9895];
    elseif alpha == 0.05
        poly_coef = [1.6597 2.4813 -16.1007 52.6998 -98.9348 91.7646 -33.3628];
    elseif alpha == 0.1
        poly_coef = [1.2917 2.4250 -14.1041 46.0326 -86.7946 80.8189 -29.4840];
    end
else
    if alpha == 0.01
        poly_coef = [2.6091 1.1854 -16.4621 63.1856 -128.0372 123.3096...
                     -45.5050 1.4378 -1.1672 -2.1843 8.4153 -9.2032...
                     3.1479 -4.7977 3.6035 -2.6765 1.0849 -0.3625...
                     12.2591 -2.5234 0.8411 0.7850 -20.5823 0.2467...
                     -0.6847 18.2815 0.6751 -6.5866];
    elseif alpha == 0.05
        poly_coef = [1.9749 1.1289 -12.2929 45.6505 -92.3587 89.5045...
                     -33.3683 1.3388 -0.8006 0.0090 0.5939 -1.0048...
                     0.2851 -4.5110 1.1262 0.9084 0.8153 -0.9854...
                     11.7294 -1.1742 -3.2329 1.7625 -18.8756 2.1281...
                     0.1723 15.5342 -0.5511 -5.2786];
    elseif alpha == 0.1
        poly_coef = [1.6552 1.2271 -11.7243 43.6253 -87.8291 84.6893...
                     -31.4176 1.2890 0.0224 -2.0585 3.2898 -2.6854...
                     0.5102 -4.8501 -0.6555 3.7550 -1.7097 0.6640...
                     14.0485 0.7875 -5.0051 1.1221 -23.9082 1.0308...
                     1.5399 20.3891 -0.5813 -7.0186];
    end   
end

%% "Prepare" the input

% Make sure Y, bounded, and bounds are column vectors
Y = Y(:);
bounded = bounded(:);
bounds = bounds(:);

% Normalize nuisance parameters to be bounded below by zero and accomodate
% lower one-sided CI
bounds(index) = 0;

Y = Y - bounds;

if type == 1
    bounded(index) = 1;
elseif type == -1
    bounded(index) = -1;
end
tool = bounded;
tool(bounded==0) = 1;

Y = tool.*Y;

Sigma = diag(tool)*Sigma*diag(tool);

% Normalize Y to have unit variances and compute correlation matrix
std_err = sqrt(diag(Sigma));

Y = inv(diag(std_err))*Y; 

Omega = inv(diag(std_err))*Sigma*inv(diag(std_err)); 

% Drop entries that are not bounded 
tool = sum(bounded(1:index-1)==0);
new_index = index-tool;

keep = (bounded~=0);
keep(index) = 1;

tool = 1:1:l;
tool = tool(logical(keep));

Y = Y(tool);
Omega = Omega(tool,tool);

% Create Y_b, Y_d, Omega_bd, and Omega_dd
l = length(Y);
ld = l-1;

Y_b = Y(new_index);

tool = ones(l)-eye(l);
tool = tool(new_index,:);

Omega_bd = Omega(new_index,logical(tool));

tool = 1:1:l;
tool = tool(tool~=new_index);

Omega_dd = Omega(tool,tool);

Y_d = Y(tool);

% Create all possible subsets
subsets = logical(dec2bin(0:2^(ld)-1) - '0');

n_ss = size(subsets,1);

%% Implement Algorithm One-Sided* or Two-Sided* using polynomial approximation
if abs(type) == 1 % One-sided CI
    s_star = 1;
    
    factor_star = 0;

    omega_star = 0;
    
    for i=2:n_ss
        Omega_bd_s = Omega_bd(subsets(i,:));
        
        tool = 1:1:ld;
        tool = tool(subsets(i,:));
        
        Omega_dd_s = Omega_dd(tool,tool);
        
        factor_temp = Omega_bd_s*inv(Omega_dd_s);
        
        if min(factor_temp) >= 0
            omega_temp = factor_temp*Omega_bd_s';
            if omega_temp > omega_star
                s_star = i;
                
                factor_star = factor_temp;
                
                omega_star = omega_temp;
            end
        end
    end
    
    if s_star > 1
        poly = omega_star.^(0:1:6);

        c = dot(poly,poly_coef);

        lb = Y_b - min(norminv(1-alpha*9/10),factor_star*Y_d(subsets(s_star,:))+c);
    else
        lb = Y_b - norminv(1-alpha);
    end

    ub = Inf;
    
    if type == -1
        ub = -lb;
        lb = -Inf;
    end
else % Two-sided CI
    s1_star = 1;
    s2_star = 1;
    
    factor1_star = 0;
    factor2_star = 0;

    omega1_star = 0;
    omega2_star = 0;
    
    for i=2:n_ss
        Omega_bd_s = Omega_bd(subsets(i,:));
        
        tool = 1:1:ld;
        tool = tool(subsets(i,:));
        
        Omega_dd_s = Omega_dd(tool,tool);
        
        factor_temp = Omega_bd_s*inv(Omega_dd_s);
        
        omega_temp = factor_temp*Omega_bd_s';
        
        if min(factor_temp) >= 0
            if omega_temp > omega1_star
                s1_star = i;
                
                factor1_star = factor_temp;
                
                omega1_star = omega_temp;
            end
        elseif max(factor_temp) <= 0
            if omega_temp > omega2_star
                s2_star = i;
                
                factor2_star = factor_temp;
                
                omega2_star = omega_temp;
            end
        end
    end
    
    if s1_star+s2_star > 2
        degree = 6;
        [i,j]= meshgrid(0:degree);
        powers= sortrows([i(:) j(:)]);
        poly_u = [];
        poly_l = [];
        for m = 1:size(powers,1)
            d = powers(m,:);
            if sum(d)<=degree
                poly_u(:, end+1) = omega1_star^d(1)*omega2_star^d(2);
                poly_l(:, end+1) = omega2_star^d(1)*omega1_star^d(2);
            end
        end
        
        c_u = dot(poly_u,poly_coef);
        c_l = dot(poly_l,poly_coef);
        
        cv_alt = norminv(1-alpha*9/20); 
        
        if s1_star > 1
            lb = Y_b - min(cv_alt,factor1_star*Y_d(subsets(s1_star,:))+c_l);
        else
            lb = Y_b - min(cv_alt,c_l);
        end
        
        if s2_star > 1 
            ub = Y_b + min(cv_alt,-factor2_star*Y_d(subsets(s2_star,:))+c_u);
        else
            ub = Y_b + min(cv_alt,c_u);
        end      
    else
        lb = Y_b - norminv(1-alpha/2);
        ub = Y_b + norminv(1-alpha/2);
    end
end

%% Return (possibly empty) confidence interval
if lb <= ub
    CI = [lb,ub]*std_err(index);
else 
    CI = [];
end

end

